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ABSTRACT 

This  paper  is  concerned  with  the  static,  one-dimensional  modelling  of  a 
semiconductor  device  (namely  the  pn- junction)  when  a  bias  is  applied.  The 
governing  equations  are  the  well  known  equations  describing  carrier  transport 
in  a  semiconductor  which  consist  of  a  system  of  five  ordinary  differential 
equations  subject  to  boundary  conditions  imposed  at  the  contacts.  Because  of 
the  different  orders  of  magnitude  of  the  solution  components  at  the 
boundaries,  we  scale  the  components  individually  and  obtain  a  singular 
perturbation  problem. 

We  analyse  the  equilibrium  case  (zero  bias  applied)  and  set  up 
approximate  models,  posed  as  singularly  perturbed  second  order  equations,  by 
neglecting  the  hole  and  electron  current  densities.  This  makes  sense  for 
small  forward  bias  and  for  reverse  bias. 

For  the  full  problems  we  prove  an  a  priori  estimate  on  the  number  of 
electron-hole  carrier  pairs  and  derive  asymptotic  expansions  (as  the 
perturbation  parameter  tends  to  zero)  by  setting  up  the  reduced  system  and  the 
boundary  layer  system.  We  prove  existence  theorems  for  both  systems  and  use 
the  asymptotic  expansion  to  solve  the  model  equations  numerically  and  analyse 
the  dependence  of  the  solutions  on  the  applied  bias.. 
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SIGNIFICANCE  AND  EXPLANATION 

This  paper  is  concerned  with  the  mathematical  modelling  of  a  diode.  This 
device  consists  of  a  semiconductor,  say  silicon,  which  is  doped  with  positive 
ions  on  one  side  and  with  negative  ions  on  the  other  aide.  We  take  the 
standard  equations  describing  hole-electron  transport  in  a  semiconductor  and 
reformulate  them  (by  scaling  of  the  dependent  and  the  independent  variables) 
as  a  singularly  perturbed  system  of  ordinary  differential  equations  (that 
means  at  least  one  derivative  of  a  dependent  variable  is  multiplied  by  a  small 
parameter,  called  perturbation  parameter)  subject  to  boundary  conditions.  The 
dependent  variables  are  the  electrostatic  potential,  the  carrier  densities  and 
the  current  densities. 

The  advantage  of  the  singular  perturbation  approach  is  that  asymptotic 
expansions  of  solutions  (as  the  perturbation  parameter  tends  to  zero)  can  be 
derived  and  used  for  the  analysis  and  numerical  solution  of  the  model 
equations. 

We  present  an  analytical  and  numerical  study  of  the  equations  (using  the 
asymptotic  expansions)  and  investigate  how  the  behavior  of  the  model  depends 
on  the  bias  applied  to  the  contacts  of  the  diode. 
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The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


A  SINGULARLY  PERTURBED  BOUNDARY  VALUE  PROBLEM 
MODELLING  A  SEMICONDUCTOR  DEVICE 

Peter  A.  Markowich*,  Chr.  A.  Ringhofer**, 
Siegfried  Selberhern***  and  Erasmus  Langer*** 


1.  Introduction 

This  paper  is  concerned  with  the  following  physical  situation.  A  semicon¬ 
ductor  (for  example  silicon)  is  doped  with  acceptor  atoms  (positive  ions)  in 
the  left  side  (p-side),  with  donor  atoms  (negative  ions)  in  the  right  side  (n- 

side)  and  a  bias  U  =  U  -  U  is  applied  to  the  contacts: 

*  AC 


cathode 
contact 

♦ 

-£  0  l  2 

_3 

The  device  is  assumed  to  have  characteristic  length  2/  =  5  x 10  cm  and  the  pn- 
junction  is  assumed  to  be  in  the  middle.  A  description  of  the  physics  of  this 
device  can  be  found  in  Ashcroft  and  Mermin  (1976)  and  R.A.  Smith  (1978). 

The  model  equations  governing  the  static  one -dimensional  case  are  (see 
de  Mari  (1968)): 


anode 

contact 


p-side 


n-side 


(a)  V  =  -(n-p-C(z)) 


(b)  .  j! '-nV  +  £.J 

n  n 

(d  p'  ^ 

n  P 


(d)  r  =  qR(n,p) 


(e)  =  -qR(n,p) 


Poisson's  equation 

electron  current  relation 

hole  current  relation 

continuity  equation  for  electrons 
continuity  equations  for  holes 
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2. 


for  -l  <  z  <  l  subject  to  the  boundary  conditions 


(a) 

+<-'>  ■  VdmAt+u* 

(anode) 

0>) 

♦  (1)  -  Vn  +  UC 

(cathode) 

(1.2) 

(c) 

n(  +  i)p(±i)  *  n* 

(d) 

n(±i)  -p(±f)  -C(  +  /)  ■=  0  , 

where  we 

require  that  ♦  ,  n  ,  p  ,  J  ,  J 

n  p 

e  Cl(t-i,i]),  *  c  C2([-i,0) U  (0,1]) 

holds 

The  set  of  full  equations  (dynamic,  two-dimensional)  can  be  found  in  Van  Roosbroeck 
(1950) . 

The  dependent  variables  (with  units)  in  (1.1)  are 


♦ 

n 


P 

J 

i 

J 


P 


potential  (V) 
electron  density  (cm”  ) 
hole  density  (an  ) 

2 

electron  current  density  (A/cm  ) 

2 

hole  current  density  (A/cm  )  . 


Of  course,  we  only  admit  solutions  fulfilling  n  >  0,  p  >  0. 

All  parameters  in  (1.1), (1.2)  (except  C(z))  and  the  temperature  (T  ~  300K) 
are  assumed  to  be  constant.  The  function  C(z),  called  doping  profile  (or  im¬ 
purity  distribution),  is  given  by 

C(z)  -  l£(z)  -h”(z)  (cm”3)  (1.3) 

where  Np  is  the  donor  density  and  N~  is  the  acceptor  density.  C(z)  is  negative 
in  [-<,0)  (p-side)  and  positive  in  (0,1]  and  is  assumed  to  jump  at  z  =  0 
(abrupt  junction).  In  this  paper  we  also  assume  that  C(z)  is  odd,  that  fc(z)| 

varies  moderately  and 


3. 


0  <  M  <  C (z)  ~ C (1)  =  1017cm-3  ,  z  e  (0,j i]  .  (1.4) 

The  physical  meaning  (and  numerical  value)  of  the  other  parameters  in  (1.1),  (1.2) 
is  given  in  Table  1. 


TABLE  1:  PARAMETERS  AT  T~300K 


Parameter 

Physical  Meaning 

Numerical  Value 

q 

elementary  charge 

10~19Asec 

e 

dielectricity  constant 

10  ^2Asec/Vcm 

^n 

electron  mobility 

103cm2/Vsec 

hole  mobility 

103cm2/Vsec 

D 

n 

electron  diffusion  constant 

or  2, 

25cm  / sec 

D 

P 

hole  diffusion  constant 

25cm2 /sec 

"i 

intrinsic  number 

1010cm"3 

D 

U  =  -^  = 
T  n 

^n 

D 

_L 

MP 

thermal  voltage 

0.025V 

Temperature-  and  space-dependent  modelling  of  these  parameters  can  be  found  in 
Langer,  Selberherr  and  Mader  (1981)  and  Selberherr  (1981) . 

R(n,p)  in  (l.l)(d),(e)  is  called  recombination  rate. 

We  take  the  Shockley-Read-Hall  (SRH)  term  (see  Langer,  Selberherr  and  Mader 
(1981)) 

2 

np"ni  -3  -1 

R(n’p)  =  r(n+p)+2rni  (cm  S6C  } 


(1.5) 


4. 


where  r = 10  ^sec  is  the  electron  (and  hole)  lifetime.  Validity  restrictions 
for  the  SRH-term  and  various  other  ways  of  modelling  the  recombination  rate 
are  given  in  the  above  cited  references  and  in  R.A.  Smith  (1978). 

The  pn-junction  is  forward  biased  if  U  =  U.  -  U  >  0  and  reverse  biased  if 

A  v 

U  <  0.  The  boundary  conditions  (1.2)  (c)  represent  thermal  equilibrium  at  z  =  i.£ 
and  (1.2)  (d)  represent  vanishing  space  charge  (the  contacts  are  assumed  to  be 
Ohmic) . 

The  system  (1.1), (1.2)  was  solved  numerically  by  Langer  (1980)  using  a 
first  order  difference  method  and  computations  for  the  corresponding  two-dimensional 
model  (and  for  MQS -transistors) and  can  be  found  in  Selberherr  (1981). 

The  approach  taken  in  this  paper  is  to  scale  (1.1),  (1.2)  in  such  a  way  that 
a  singular  perturbation  problem  is  obtained.  This  is  suited  to  the  problem  since 
the  solutions  exhibit  internal  layer  structure  (at  z=0)  due  to  the  jump-discon¬ 
tinuity  of  the  doping  profile  (and  since  the  coefficients  in  (1.1)  are  of  different 
orders  of  magnitude  (see  Table  1  and  (1.4)).  Moreover  the  oddness  of  C(z)  will 
be  used  in  order  to  obtain  an  equivalent  problem  on  the  half-interval  [0,2] 

(which  is  of  boundary  layer  type) . 

The  advantage  of  this  procedure  is  two-fold.  Firstly,  asymptotic  expansions 
of  solutions  (as  the  perturbation  parameter  converges  to  zero)  can  be  obtained  and 
secondly  the  problem  can  be  solved  numerically  using  already  developed  (high  order) 
methods  for  boundary  layer  problems  (for  example,  polynomial  collocation  as  used  in 
this  paper).  Vasileva  and  Butuzow  (1978),  Vasileva  and  Stelmakh  (1977)  and  D.  Smith 
(1980)  proceeded  in  this  way  obtaining  asymptotic  expansions  for  a  simplified  problem 
(they  assumed  that  the  recombination  rate  R  is  a  known  function  only  depending  on  z 

The  paper  is  organized  as  follows.  In  Section  2  we  perform  the  scaling  and  re¬ 
formulate  (1.1),  (1.2)  as  boundary  layer  problems  on  the  half-interval.  Section  3 


5. 


Is  concerned  with  the  case  U=0  and  with  zero-current -density  models.  In 
Section  4,  we  derive  asymptotic  expansions,  demonstrate  some  qualitative  prop¬ 
erties  of  the  solution,  and  explain  the  numerical  results.  In  Appendix  A  we 
give  proofs  of  Theorems  on  the  zero-current-density  models  and  Appendix  B  is 
concerned  with  collocation  methods  for  singular  perturbation  problems. 


and  get  the  system 
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*"  -  £-"^^(n  -p  -D (x) ) 


n  \|r'  +  J 

a  s  tD  n 
n  s 


p«  .  „  p  \|»*  _  ‘iL-  j 

*s  »  8  rD  p 

P 


-1  <  x  <  1 


(2.4) 


ul  2 

i  W^cTIT^ 

4  n± 

ns+ps+2C(i) 


J'  -  -J' 
Po  “ 


with  boundary  conditions 


ni 

<7T777> 


(b)  +s(l) 


°s(1)  UC 

in-\-  +  ”7 


(2.5) 


(c)  ns(±l)p8(±l)  «  (-L)2 


(d)  ns(±l) -P8(±l)  -  (±D  -  0 


We  set 


„  U_e 

,  2  T 


i2C(i)q 


2  _  V  « 

Y  IT  *• 


(2.6) 


Using  the  values  for  the  parameters  given  in  Section  1  and  dropping  the  sub¬ 


script  s,  we  obtain 


7  . 


(a) 

(b) 

(c) 

(d) 

(e) 


aV  =  n-p-D(x) 


n '  =  n'l'1  +  J 


p’  =  -  p'K  -  J 

P 

T.  _  I  np-y^A4 
"  ■  4 

j<  =  1  -hP-V-V 

P  4  n+p+2y2A2 


for  -1  <  x  <  1  subject  to  the  boundary  conditions 


*(-l)  =  , 

= 

A2 

UA 

(a) 

in 

■4  ml' - 

P(-1) 

+  — 
UT 

(b) 

*(1)  =  *+  , 

V  * 

in 

aill. 

yV 

(c) 

n(il)p(±l)  = 

yV 

(d) 

n(  +  l)  -p(±l) 

-  (±D  = 

=  0 

(2.6) 


(2.7) 


?  -6  2 

where  A  =  0.4x10  ,  y  =0.25.  All  quantities  in  (2.6),  (2.7)  are  dimensionless. 

The  problem  (2.5), (2.6)  will  be  regarded  as  singular  perturbation  problem 

with  perturbation  parameter  A.  Physically,  decreasing  A  (A  >  0+)  means  that 

2 

the  concentration  of  impurities  increases  (C(i) ->-<»).  Note  that  y  is  independent 
of  C(i). 

So  far,  we  have  not  heavily  used  the  oddness  of  C(z)  (except  for  (2.1)  and 
(2.7)  (d)).  For  the  performed  scaling  we  only  have  to  require  that  |c(z)|  is  of 
the  same  order  of  magnitude  for  all  z  e  [-i,f].  Now  we  will  use  the  oddness 
assumption  in  order  to  reduce  the  problem  (2.6),  (2,7)  to  the  half-interval  [0,1], 
thus  making  the  Internal  layer  at  x=0  (called  depletion-layer  or  space-charge 
region  in  this  context)  to  a  boundary  layer. 


8. 


Since  only  'K'  ,  enter  the  equations  (2.6),  the  substitution  ^  (x)  =  ^(x)  - 

* 

- - -  does  not  change  (2.6).  Boundary  conditions  for  are 

\|/  _\|* 

V*(l)  ='!»*  =  ~~2~  »  'M“1)  =  (2-8) 

and  >[/  (-l)-'l'  (l)='tr  -'k  holds.  Therefore  we  can  assume  (without  loss  of  generality) 
that 

(uc=-uA)  (2-9) 

T 

holds.  Moreover,  given  a  solution  ('I'.n, p, Jfl  ,  J  )  of  (2.6),  (2.7)  we  immediately 
verify  that 


't'j(x)  =  -4(-x) 

nx(l)  =  p(-x) 

Pt(x)  =  n(-x)  (2.9) 

J  (x)  =  J(-x) 
nl  P 

J  (x)  =  J(-x) 

P1  n 

constitutes  another  solution,  since  D  is  odd  and  since  the  SRH  recombination 
rate  (1.5)  is  symmetric  in  n  and  p.  Therefore  we  investigate  (2.6)  on  [0,1] 
subject  to  the  boundary  conditions 

(a)  *(1)  =  *+  ,  *(0)  =  0 

(b)  n(l)p(l)  =  yV  ,  n(l)  -p(l)  =  1 

(2.10) 

(c)  n(0)  =  p(0) 

(d)  J  (0)  =  J  (0) 

n  p 

and  continue  the  solution  to  [-1,1]  according  to  (2.9). 
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We  have  to  expect  a  boundary  layer  (as  A  >  0  +)  at  x  -  0  in  the  solution 
components  ^  ,  'V'  ,  n,  p  (which  are  the  fast  variables),  the  slow  (unperturbed) 
components  will  be  uniformly  smooth  (in  A),  no  boundary  layer  will  occur 

at  x=l  since  the  reduced  solution-manifold  (A=0),  which  is  given  by 
n-p-D(x)  =  0,  satisfies  the  (reduced)  boundary  condition  (2.10)  (b),  given  by 
p(l)  =  0,  n(l)  =  l.  A  proof  for  this  will  be  given  in  Section  4. 


3.  Zero-Current  Approximations 

We  are  now  concerned  with  the  calculation  of  the  built-in-potential  (no 

voltage  is  applied,  U  = U  =  0) .  For  simplicity  we  assume  D  =  1  in  this  Section. 

A  L/ 

Setting  J  sJ  =0  we  get  from  (2.6)(b),(c) 
n  p 


n  =  Ae^  ,  p  =  Be  ^  ,  A,B  e  R  (3.1) 

and  from  (2.6)(d),(e)  we  obtain  AB=y^A^.  The  boundary  conditions  (2.10(a),  (b),  (c) 
imply  A  =  B=y'7^  and 


'Kl)  =  in 


1  4-  \  1  +  4-y‘Ss^ 

2yh2 


:  *°(A) 


\|/(0)  =  0 


(3.2) 


Inserting  (3.1)  into  (2.6) (a)  gives  the  equation  for  the  built-in-potential 


2  2.  2  ^  -\!/ 

Ay"  =  yV(e  -e  v  )  -  1  ,  0  <  x  <  1 


(3.3) 


subject  to  the  boundary  conditions  (3.2). 

This  problem  is  slightly  nonstandard  since  VjJ(A)  ~  in  -j ^  >  «  as  A  >  0  + 

Y  A 

and  since  the  reduced  equation  (3.3)  (A  =  0)  is  inconsistent. 

We  obtain  the  following  existence  result: 


Theorem  3.1.  For  A  sufficiently  small,  the  problem  (3.2),  (3.3)  has  a  locally 
unique  solution  4^  fulfilling 


10. 


_  c 

^(x.X)  =  Vj (A)  (1-^(2)  +  0 (e  *)) 

where  C  >  0,  V  (0)  ■  1,  V  (y)  t  0  as  y  •>  »  and 

A  A 

0  <  VA(y)  <  exp((-l  +  6)y)  ,  y  >  const,  in  in 

holds  for  some  0  <  6  <  1. 

The  proof  is  given  in  Appendix  A. 

The  corresponding  carrier  densitits  n,p  are  given  by  (3.1): 

nQ(x,A)  =  exp (-  ^(MV^Hl+OfrV)) 

P0(x,7O  =  vVexp(  ^)V^(|))(1+0(yV))  . 


(3.4) 

(3.5) 


(3.6)  (a) 

(3.6)  (b) 


Since  U*  2Ut('I'^(A)-'I'+)  the  pn- junction  is  forward  biased  for  <  ty+(A)  and 

reverse  biased  for  >  ^(A) . 

Theorem  3.1  encourages  us  to  consider  approximate  problems  for  (2.6),  (2.10) 
(with  ~  Vj|(X))  obtained  by  setting  =  ®  in  (2.6)(b),(c)  and  ignoring 

(2.6)(d),(e)  (equilibrium  models) .  For  the  two-dimensional  case,  this  was 
suggested  by  Polak  et  al  (1981) . 

One  of  the  boundary  conditions  in  (2.10) (b), (c)  has  to  be  ignored 
(since  (2.10)  (b),  (c)  are  consistent  with  J^sJ^sO  iff  ^  =  'k+(A)) . 

At  first  we  drop  the  condition  n(l)p(l)  =  yV  (thermal  equilibrium  at  the 
contacts)  and  get  from  (3.1), (2.10) (b), (c) 


'I'(O)  =  0  ,  *(1)  - 


(3.8) 
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For  =  (3.7),  (3.8)  gives  (3.3),  (3.2).  Because  we  require  n  >  0,  p  ^  0, 

the  model  (3.7),  (3.8)  is  relevant  only  for  >  0. 

We  get  as  a  generalization  of  Theorem  3.1: 

Theorem  3.2.  Assume  that  A  is  sufficiently  small  and  that 

0  <  C.  <  *.  <  — l-- -  (3.9) 

2A 

holds  for  f  <  cj  <  1,  0  <  &  <  1  independent  of  A  ,  .  Then  the  problem  (3.7), 

(3.8)  has  a  locally  unique  solution  =  'l' (x, A,^  )  fulfilling 

e  + 

=  *+(l-V+  (~)  +  0(exp(-^~%)  ,  A  *  0  +  (3.10) 

uniformly  in  where  V,  (0)  =  1,  V  (y)  i  0  as  y  •>  <*■,  and 

0  <  (y)  <  exp ((-1  +  £)y)  ,  y  >  const.  (3.11) 


holds . 

The  proof  is  given  in  Appendix  A. 

The  carrier  densities  ng  ,  p&  are  calculated  from  (3.1),  (3.7)  giving 


♦ 

e 

_  e _ 

e  2sinh(  'l'+)  * 


2sinh(  ) 


(3.12) 


For  the  existence  proof  (and  for  n,p  >0)  without  uniform  asymptotics  (in  V  ) 
>  0  is  sufficient  . 

The  width  d  of  the  depletion  layer  of  ^  is  estimated  using  (3.10): 

d  <  const.  A  in |-^-|  as  A  +  0+  (3.13) 

uniformly  for  'l'+  satisfying  (3.9)  where  const.  ~  1. 


The  width  of  the  layer  depends  logarithmically  on  ♦  for  fixed 

Figures  1  and  2  show  the  computed  potential  lij.'l'  and  the  hole  density 
respectively  (on  the  interval  [0,0.0012])  for  =  10,  which  corresponds  to 

U  ~  0.306V  (forward)  bias  since  (for  the  parameters  given  in  Section  1) 

♦?0O  »  16.118. 

~T 


The  second  equilibrium  model  we  investigate  is  obtained  by  dropping  the 
symmetry  condition  n(0)  =  p(0).  From  (2.10)  (b)  we  get 

n(l,M  =  |(1  Wl  +  4yV  )  ,  p(l,M  =  j(-l  +n/i  +  4YV)  (3.14) 

-  \Jr 

T-f  + 

and  therefore  A  =  n(l,70e  ,  B=p(l,>)e  .  Poisson's  equation  then  reads 

(♦-♦  )  -  (♦-'If.) 

=  n(l,A)e  -p(l,-\)e  -1,  0<x<l  (3.15) 

subject  to  the  boundary  conditions  (3.8). 

We  get 

Theorem  3.3.  For  every  ♦  e  B  and  ~h  <  C0|f+)  the  problem  (3.15),  (3.8)  has 
a  locally  unique  solution  ♦  fulfilling 

*a  =  *+(l  -W^  (^)  +  0(yV))  ,  *  *  0+  (3.16) 

"T* 

where  W^.  (0)  =  1,  (y)  10  as  y  >  °°  and 

0  <  W.  (y)  <  const.  exp((-l+6)y)  ,  y  >  0  (3.17) 

v+ 

holds.  (For  the  proof  see  Appendix  A.) 

Figures  3  and  4  show  the  computed  curves  U  ♦  and  p  respectively  (on 
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In  order  to  set  up  the  third  possible  equilibrium  model  the  vanishing-space- 
charge  condition  n(l)  -p(l)  =  l  has  to  be  neglected.  Since  the  analysis  for  this 


is  quite  similar  to  the  analysis  of  (3.7),  (3.8),  we  do  not  give  details. 

A  comparison  of  numerical  results  obtained  for  the  equilibrium  models  and  the 

full  nonequilibrium  model  (see  next  section)  shows  that  the  agreement  is  excellent 

-9 

as  long  as  U  <  0.3V  (compare  Figures  1  and  15)  since  then  J=J  +J  <5x10 

n  p  ~ 

holds.  The  agreement  gets  rapidly  worse  as  U  increases,  since  J  increases  ex¬ 
ponentially  with  U  (compare  Figures  3  and  8) . 

4.  The  Nonequilibrium  Problem 

At  first  we  prove  an  a  priori  estimate  on  the  number  of  electron-hole  (carrier) 
pairs. 

Theorem  4.1.  Assume  that  the  problem  (2.6)  and  (2.10)  has  a  solution  for  which 

n  >  0,  P  >  0  holds  and  Jn  ,  do  not  change  sign  in  [0,1] .  Then  Jn  >  0, 

J  >0,  J  ^0,  J  #0  holds  for  forward  bias  (*!'  <  ^(A))  and  J  <  0,  J  <0, 
p  —  ’  n  '  p  +  +  n~  p  — 

J  #  0,  J  #0  holds  for  reverse  bias  ('I'  >  V?(A)).  Moreover,  the  estimate 
n  p  ■+■  + 

Y^A^  exp(  -  -|— ^-)  <  n(x)p(x)  <  ySv^exp  (-^~4  >0<x<l  (4.1) 

holds . 

The  assumption* of  the  Theorem  are  physically  reasonable  since  n,p  are  the 
(scaled)  numbers  of  carriers  and  since  the  currents  do  not  change  direction  in  the 
device. 

Proof.  It  is  more  convenient  to  work  with  the  problem  (2.6),  (2.7)  (with 


We  get  from  (2.6)(b),(c) 


n  *=  (E  +  T  e“  (s)ds)e  * 

*■^-1  n 

(4.2) (a) 

p  ■  (f  -  . 

_l  P 

(4.2)  (b) 

-  v(/  \jr 

Obviously  n(-l)*Ee  +,  p(-l)«Fe  +  holds. 

From  (2.7)  (c)  we  conclude  and  (2.7)  (d)  gives 

-  'Sf  \lf 

+  + 

Ee  -  Fe  =  -1. 

Prom 

these  equations  and  from  n  >  0,  p  >  0  we  calculate 

E  =  |e  %-l  Wi  +4yVS 

(4.3)  (a) 

F  -  j  e  *(1  +A  +  4y\4  ’)  . 

(4.3)  (b) 

We  define 

X  -  J  «*  ^(8)Jn(s)ds 

(4.4) (a) 

Y  -  f  e  '^J  (s)ds 

J-1  p 

(4.4)  (b) 

and  get  from  the  boundary  conditions  (2.7(c),(d)  at  x=l: 

(E  +X)  (F  -  Y)  =  A* 

(4.5) (a) 

'if  -  'if 

e  +(E  +X)  -  e  +(F-Y)  =  1  . 

(4.5) (b) 

This 

system  of  equations  has  the  solutions  X  * Y  *  0  and 

X-F-E,  Y  =  X  . 

(4.6) 

The  trivial  solution  X«Y“0  Implies  J^J^sO  (since  Jn  *  are  not 
allowed  to  change  sign)  and  this  Implies  np=  A4  and  't'  =  ^  (A )  .  Therefore 

”T  T 

we  take  (4.6).  X* Y  Is  also  an  immediate  consequence  of  the  syrmnetry  (2.9). 


1 1 
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A  simple  calculation  shows  that  X  >  0,  Y  >  0  iff  <  4^  (A)  (forward  bias) 

and  X  <  0,  Y  <  0  iff  'k  >  4^  (A)  (reverse  bias).  This  proves  the  sign-statement 

on  J  ,  J  . 
n  p 

Also 

n(x)p(x)  =  (E  +  f  e  (s)ds)  (F  -  T  e  (s)ds)  (4.7) 

J_1  n  p 

holds .  Therefore  we  get  for  forward  bias 

E2  =  E  (F-X)  <  n(x)p(x)  <  (E  +  Y)F  =  F2  (4.8)  (a) 

and  for  reverse  bias 

F2  =  (E  +  X)F  <  n(x)p(x)  <  E(F-Y)  =  E2  .  (4.8)  (b) 

Using  (4.3)  and  (3.2)  gives  (4.1). 

Expressing  (4.1)  in  the  unsealed  variables  n,p  (2.2)  (a)  gives 

n? exp (  --^-4  <  n(z)p (z)  <  n2exp(-^)  ,  -i  <  z  <  i  .  (4.9) 

This  estimate  was  anticipated  by  de  Mari  (1968).  There  is  numerical  evidence  that 
for  forward  bias  n(z)p(z)  is  close  to  the  upper  bound  while  for  reverse  bias  it 
is  close  to  the  lower  bound. 

We  remark  that  the  explicit  form  of  the  recombination  rate  R(n,p)  was  not 
used  for  the  proof  of  Theorem  4.1  (except  for  excluding  the  equilibrium  case 


Now  we  assume  that  D  e  (^([0,1])  and  that  the  solutions  of  (2.6),  (2.7)  have 


an  asymptotic  expansion  given  by 


+  (*.*) 
n(x,7v) 
p(x,A) 


Jn(x,A) 


Jp (x,A) 


^  (x) 

tJj  (x) 

Pj  (x) 

?  (x) 

"j 

Vx> 


*,(r) 
tlj  (t) 
pj(t) 

Jn.(T) 

J 

J  (r) 
pj 


V” 

Sj(«) 

PjO) 

J  (a) 

"j 

J  (cr) 

Pj 


r  =  ->a  =  — 


(A. 10) 


where  ,  n  ,  p  ,  J  ,  J  are  smooth  on  [0,1],  'i'.(r),  n  (t),  p,(t), 
j  j  J  nj  Pj  J  J  J 

J  (t),  J  (r)  and  $.(cr),  n.(cr),  p.(o),  J  (cr),  J  (cr)  are  defined  on 
°j  Pj  j  j  J  nj  Pj 

[0  ,«>)  and  decay  (exponentially)  to  zero  as  r  >  »  and  a  >  «>  respectively. 
The  functions  marked  with  are  the  boundary  layer  terms  decaying  from 

the  left  boundary  and  the  functions  marked  with  '  A  '  are  the  boundary  layer 
terms  decaying  from  the  right  boundary. 

Inserting  (4.10)  into  (2.6),  comparing  coefficient,  setting  ~K~  0  and 
using  that  the  boundary  layer  terms  vanish  asymptotically  away  from  the 
corresponding  boundary  gives  the  reduced  problem 


(a) 

0  = 

n0-p0-D(x) 

(b) 

n0  = 

=  “0*0  +  5n0 

(c) 

Po  = 

(d) 

J' 

n_ 

i  Vo 

4  ~  ~ 

0 

”o  +  po 

(e) 

J' 

po 

1  Vo 

"  4  ~  ,  ~ 
n0  +  p0 

0  <  x  <  1 


(4.11) 


the  (zeroth  order)  left  boundary  layer  system 


*0  •  no  “  po 

V  <W0))^o 


Po  =  “  (Po  +  Po(0))'l'o 


J  *  0 
n0 


J  =0 
p0 


0  <  T  <  * 


(4.12) 


and  the  fceroth  ordei)  right  boundary  layer  system 


(a) 

*  A 

*0  "  n0 

• 

A  A 

(b) 

n0=  (ni 

(c) 

P0  =  -  (i 

• 

A 

(d) 

J  *  0 

no 

• 

A 

(e) 

J  =0 
Pn 

0  “O' 

A  iv 

»o+,,o' 


0  <  a  <  °° 


(4.13) 


where  1  .  ’  denotes  differentiation  with  respect  to  r  or  o  respectively. 
The  boundary  conditions  (2.10)  give 


V0)  +  V0)  =  0 

nQ(0)  +  nQ(0)  =  pQ(0)  +  ?0(0) 


Jn  <°>  +  Jn  (°>  =  Jn  (0)  +  J  (0) 

n0  no  P0  P0 


(4.14) 


^0(1)  +  MrQ(0)  =  *+ 


n0(l)  +  nQ(0)  =  1 
P0(D  +  P0(0)  =  0 


(4.15) 


}  o. 
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since  (2.10) (b)  reduces  to  n(l) =  1,  p(l)  =  0  for  \  =  0.  The  systems 
(4.11),  (4.12),  (4.13)  were  obtained  by  proceeding  similarly  to  O’Malley  (1978), 
Vasileva  and  Butuzov  (1978)  and  Vasileva  and  Stelmakh  (1977). 

From  the  decay  of  the  layer  terms  and  from  (4.12) (d), (e),  (4.13) (d), (e), 
we  immediately  get 


(4.16) 


Jq  ,  .1^  are  the  slow  (unperturbed)  components  which  do  not  exhibit  layer  behaviour. 
Now  we  show  that  no  boundary  layers  at  x=  1  occur. 


Integrating  (4.13)  (b),  (c)  gives 

nQ(o)  =  S0(l)(e  °  -1),  P0(<O  -  P0(D(e  0  -1)  (4.17) 


v°>  ~  ■  V0) 

From  (4.15)  (b)  we  conclude  ng(l)e  =1  and  Pg(l)e  =0  and  therefore 

(1) =  0.  This  gives 


P0(c)  h  0  .  (4.18) 


(4.11)(a)  gives  Hq(1)  =  1  and  therefore  'I'q(O)  =  0  holds.  From  (4.13)  (a)  and 

A 

from  the  decay  of  (a)  we  conclude 


$Q(a)  a  0  ,  nQ(a)  h  0  .  (4.19) 


Boundary  layers  of  order  zero  can  only  occur  at  x  =  0. 

Integrating  (4. 12) (b), (c)  gives 

_  ^(r)  _  'Mr) 

nQ(r)  =  nQ(0)(e  -1)  ,  pQ(r)  ■  pQ(0)(e  -1)  .  (4.20) 


From  (4.14)  (b)  we  conclude 
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v°>  ~  -  V°> 

nQ<0)e  =  pQ(0)e 


(4.21) 


and  since  'l'g(O)  = -'l'g(O)  ((4.14)  (a))  we  find  the  missing  boundary  condition 


for  the  reduced  system: 


tig  (0)  e 


-  v°> 


P0(0)e 


V°> 


(4.22) 


The  reduced  system  of  equations  can  now  be  obtained  in  closed  form  from  (4.11)  by 
substituting  pQ=nQ-D(x)  into  (4 . 11)  (b) - (e)  and  by  eliminating  from 

(4.11)  (b),  (c) .  We  drop  the  subscript  0  and  get 


(a) 


(b) 


(c) 


(d) 


= 


n(D 1  (x)  +  J  -J  )-D(x)J 
_ n  p _ n 

2n-D (x) 

D '  (x)  -  (J  +J  ) 


2n-D  (x) 


Jn  "  4 


=  i-  n^n-D^xjj, 
2n-D (x) 


>  0  <  x  <  1 


(4-23) 


j.  -  In(n-D(x)) 
p  2n-D (x) 


(e) 

p  =  n  -  D  (x) 

the  boundary  conditions 

(a) 

y-N 

o 

V-/ 

II 

o 

(b) 

(n(O)-D (0) )e2  *(0) 

(c) 

II 

✓“N 

H 

(d) 

n(l)  =  1  . 

=  n(0) 


(4.24) 


From  (i . 1 2) (a) ,  (4.20)  we  get  the  boundary  layer  equation 


+  =  n(0)e  *  -  (n(0) -D (0) )e”^-D (0)  ,  0  <  x  <  -  (4.25; 

f(0)  =  -5(0)  ,  50*)  =  0  .  (4.26) 

Higher  order  terms  of  the  expansion  (4.10)  can  be  obtained  analogously,  however, 
we  will  not  do  so. 

At  first  we  prove 

Theorem  4.2.  Assume  that  D(0)  >  0,  p(0)  =n(0)  -D(0)  >  0.  Then  the  boundary  layer 
equation  (4.25),  (4.26)  has  a  unique  monotonely  increasing  solution  9(t),  ful¬ 
filling 

-C&exp(-' J  (n(0) +p(0))(l-8)r)  <  ^(r)  <  -C"lexp (- n/  (n(0)  +  p(0))  (1 +B)r)  (4.27) 

on  [0,a=]  for  any  6  >  0  with  CL  >  0. 

o 

Proof .  The  substitution  <p=  gives  the  problem 

cp  =  p(0)e  ^  -  n(0)e_  ^  +  D(0)  <.4.28)(a) 

9(0)  =  5(0)  ,  <p(°o)  =  0  .  (4.28)  (b) 

For  D(0)  >  0,  n(0)  -D(0)  >  0  we  get  5(0)  >  0  from  (4.24)  (b).  Applying  Fife's 
(1973)  Lenina  2.1  to  (4. 28) (a) ,  (b)  yields  the  result. 

If  D(0)  *  0  holds,  a  locally  unique  solution  of  (4.25), (4.26)  is  given  by 
^  =  0  which  implies  n;p=0. 

No  boundary  layers  occur  if  the  junction  is  not  abrupt  (see  also  D.  Smith 
(1980)  for  a  simplified  model). 


The  width  d  of  the  boundary  layer  (depletion  layer)  of  '<  ran  be  calculated  by 
taking  C  =  'I'(O)  and  by  expressing  n(0),  p(0)  in  terms  of  'V(O). 


d  ~  A 


'  2  '?’(°)+i ' 

D(0)  2 - — 

\  e2^0)-!/ 


-h 


£n 


e2*<°hl 


,2^(0).! 


%  \ 


(4.29) 


(4.29)  should  be  compared  to  (3.13)  taking  into  account  the  Theorem  4.3. 


Assuming  the  validity  of  the  expansion  (4.10),  p(0)  >  0  (which  implies 

n(0)  ^D(O)  n  0)  is  necessary  for  p,n  ^  0  if  A  is  sufficiently  small  since 


n(0,A)  =  n(0)  +  n(0)  +  0(A)  =  n(0)e“  3?(0)  +  0(A) 


P  (0,A)  =  p  (0)  +  9(0)  +  0(A)  =  p(0)e 


*<°>  +  0 


(A) 


(A. 30)  (a) 
(4.30) (b) 


holds.  p(0)  =  0  is  impossible  because  of  (4.24)(b). 

Since  4'+  =  'l'+(A)  -~j-holds,  the  reduced  problem  (4 . 23) ,  (4 . 24)  can  be  regarded 
as  dependent  on  A  for  U  <  Uq  .  We  now  prc»  »  an  existence  result  for  this  case. 

Theorem  4.3.  Assume  that  D  e  C^([0, 1])  and  D(x)  >  M  ^  0  on  [0,1|  and  that 

U<UQ,  UQ  <  2*°(A)Ut  (4.31) 


holds .  Moreover  if 


£*0 


(4.24) 


holds,  then  the  reduced  problem  (4.23),  (4.24)  has  a  locally  unique  solution 
fulfilling 


26. 


t 


(a) 

5(x)  «  ♦ .  +  in  D  (x)  +  0(-y) 

>  00 

(b) 

n(x)  =  D  (x)  +  0(“7) 

*+ 

'jf  .  ■>  °° 

T 

(c) 

p(x)  =  0(-^) 

*  “ 

(4.25) 

(d) 

Jn(x)  -  0(i) 
y+ 

>j/+  >  00 

(e) 

J  (x)  =  0(-4) 
p  ♦; 

>  CO  . 

Proof.  We  substitute  <P=  •  >  9e*-  ^4.  ~  ip  an<^ 

-  \|/  y+  + 

get  from  (4.23) (b), 

(4.24)  (b), (c) 

D '  -  ( J  +J  ) 

~  ~  .2  n  p 

<p'  ■  -0p  +  <P.)  z 

2n-D 

(4.26) (a) 

n(0)exp(-“ - )  -  (n(O)-D(O))  * 

9(0)-Hp+ 

0 

(4.26) (b) 

5(1)  =  0  . 

(4.26)  (c) 

We  rewrite  the  reduced  problem  in  operator  form  F  6)  =  0  where 
F  :  0  c  (C1  CIO,  11))^  *  (C((0,1]))4  x  R4,  0  open  and  z  =  (v.n.J,  ,  Jp) : 


F  (z)  -  (-p 
<P. 


.2d'-(W  °(II'4J.-Jp)-pJn  J,  .la&a 

+  (*  +  ♦+>  2-_D  .  -  2~.D  •  »  4  2i-V 


j,  +  I  5&21  ,  5(i,  ,  5(i).  l,  j(0)-Jn(0), 

P  **  ^ 


(4.27) 


2n-D 


n(0)exp(-^-^ - )  -  (n(O)-D(O)))  . 

<p(0)+<p+ 


The  spaces  are  equipped  with  their  natural  norms  (see  Appendix  A) . 


27. 


We  define 


z0  =  (O.D.O.O)  e  ft 


(4.28) 


and  evaluate  F^  at  z^: 


F  (O  *  (-—-,  0, 0,0,0, 0,0, D(0)exp( )) 

tD  U  U 


■<p+v  0 


(4.29) 


This  yields  consistency 


2 

||F  (z  )||  <  const.  <p+  as  tp+  >  0  + 

<p+  0 


(4.30) 


We  calculate  the  linearization  F'  (zQ)  and  investigate  the  equation 

M- 

F'  (z  )  =  (f  o>)  f  e  <C([0,ll))*,  co  e  H4.  So  we  obtain  the  linear  bound  ary -value 
<p+  0 

problem: 

r2  2 

D'  2 

-2cp+  -5-  2<p+  ^2  D  D 


0  -1 


0  0 


y  +  f,  0  <  x  <  1 


(4.31) (a) 


28. 


yl  (!)  =  ^ 

y2(l)  =  u>2 
y3 "  y4  (0^  =  “3 

i-2|2I  exp(-|-)yi(0)  +  (exp(-^-) -l)y2(0)  =  u>4  . 

«p+  +  + 

The  perturbed  problem 


(4.31)  (b) 


29. 


vj  =  -(in  DJ’Vj,  -v3 


4  1 


■  A. 

V  S  a  “•  V 

3  4  1 


v1(l)  =  0  ,  v^O)  =  0  ,  v2(0)  =  v3(0) 


(4.34) (a) 


(4.34) (b) 


has  only  the  trivial  solution  v  =  0.  This  holds  Iff 

y"  +  (in  D)'y'-|-y  =  0  ,  y'  (0)  =  y' (1)  =  0 


(4.34)  (c) 


has  only  the  trivial  solution  y  =  v3sO,  which  follows  from  the  maximum  principle 
since  <  (in  D) '  <  M2  holds. 

The  contraction  mapping  theorem  assures  the  unique  solvability  of  (4.31)  (for 
<p+  sufficiently  small)  for  all  f  e  (C([0, 1]))\  to  e  and  the  stability 
estimate 


Ibll 


to.i] 


I[0,1J  -  Kl' 


![0,1] 


4 

z 

i=l 


1^1) 


(4.35) 


where  is  independent  of  <p+  . 

An  easy  calculation  shows  the  uniform  Lipschitz  continuity  of  locally 

around  z. 


11*1  <5)-F'  (v)||  <K  ||u-v||  (4.36) 

with  ||u-Zq||,  ||v-Zq]|  <  p  and  K2  is  independent  of  <p+  .  Applying  Keller '  s  (1975) 

theory  we  conclude  that  the  equation  F  (z)  =  0  has  a  locally  unique  solution 

v+ 

z=z(«p+)  if  <p+  <  c(K^  ,  K2)  which  fulfills 


30. 


|*(9+)-*0ll[0  I] 


r> 

+  ||z' (<P+)-z'Hro  <  const. 


(4.37) 


This  translates  immediately  into  the  existence  statement  given  In  Theorem  4.3 
1 


(using  <p+  =  “o 


*;w-2s 


).  (4.25(b)-(e)  follow  from  (4.37)  and  (4.25) (a)  follows 


from  (4.23) (a). 


If  Dil  then  the  first  component  In  F^  (z^)  is  zero  and  instead  of  (4.30) 


we  get 


II*  =  exp("m_)  ’  (4.38) 

Therefore  the  0 (-4:) -terms  in  (4.25)  can  be  substituted  by  0(exp(-2  'if  )) 

and  the  bound  in  (4.24)  can  be  chosen  larger. 

Questions  that  arise  are  whether  the  asymptotic  expansion  (4.10)  actually 
represents  a  solution  of  (2.6), (2.10)  and  whether  (4.10)  is  uniform  in  4»+  in 
some  sense  (this  is  particularly  interesting  in  context  with  Theorem  4.3).  For 
simplified  pn- junction  models  the  first  question  was  positively  answered  by  D.  Smith 
(1980),  Vasileva  and  Butuzow  (1978)  and  Vasileva  and  Stelmakh  (1977).  For  a 
validity  proof  of  the  asymptotic  expansion  (4.10)  (for  ^+  =  0(1)  as  A>0+)  one 
has  to  proceed  similarly  to  the  references  cited  above  (particularly  Vasileva 
and  Butuzow  (1978)). 

Our  numerical  experiments  show  that  the  expansions  are  valid  in  the  sense 
that  a  numerical  solution  of  (2.6),  and  (2.10)  can  be  obtained  by  using  the  sum 
of  the  reduced  solution  and  the  boundary  layer  solution  as  initial  guess  for  the 
numerical  scheme  for  (2.6)  and  (2.10)  and  this  obtained  solution  is  close  to  the 
initial  guess. 

For  the  following  computations  the  parameter  values  from  Section  1  and 


D(x)sl  were  taken. 


31. 


Figures  5  through  7  show  the  reduced  solutions  ,  n  ,  Jn  ,  for  ^=0 

which  correspond  to  a  forward  bias  U=  0.806V.  The  corresponding  full  solutions 

U  '!'  ,  n,  p  (of  (2.6),  (2.10))  are  shown  in  Figures  8  and  9.  J  and  J  are  not 
i  n  p 

shown  since  they  differ  inappreciably  from  .  Figures  10  and  11  show 

Jn  ,  Jp  for  U=  1.506V  (^=-14)  and  for  U=  0.306V  (v|f+=10)  respectively. 

We  get  the  following  Table  2  for  the  total  current  J  =  J  +J  ~J  +J 

n  p  n  p 


U 

J 

10 

0.306V 

4,4  x  10-9 

0 

0.806V 

0.92 

-14 

1.506V 

108.2 

Table  2:  Current -Voltage  Relation 


The  characteristic  curve  J  =  J (~~)  behaves  exponentially  (see  A.  Smith  (1978), 

UT 

Section  7.9). 

The  computations  demonstrate  that  $(0)  0  +  as  ^  >-»  and  therefore 

the  jump  of  the  boundary  layer  of  'I'  tends  to  zero.  Since 

e  *(0) 

n(0)  =  -D(0)  — - -  (4.39)  (a) 


p(0)  -  D(0)  — -1 -  (4.39)  (b) 

e 

hold,  the  layer  jumps  of  n  and  p  tend  to  2^21.  .  This  is  illuatrated 

by  the  Figures  12  and  13  which  show  the  solutions  and  n,  p  (of  (2.6),  (2. 10) ) 

respectively,  for  '|>+*  -14  (U  *  1.506V).  No  layer  in  is  visible  and  the  layers 


isd-in 


* 


j 
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of  n,  p  are  of  approximate  height  0.5  (see  also  Figure  14  which  shows  the 
corresponding  boundary  layer  solutions  n  ,  p) . 

Figures  15  and  16  show  the  solutions  and  n  respectively  with 

U  =  0.306V (♦  =  10) .  The  corresponding  hole  density  p  is  nearly  constant  zero 
(indistinguishable  from  the  x-axis).  'i'  and  n  are  (disregarding  the  layers) 
nearly  constant.  The  explanation  for  this  behaviour  is  given  by  Theorem  4.3. 

Since  D(x)=l  the  estimate  (4.25)  can  be  strengthened  (see  remark  to  Theorem  4.3): 


iiH-i 


[0,1] 


+  ln-1 


[0,1] 


+  !|p|l[0>1]  +^J[0,1]  +^[0,1]  -  const- (e*P(-2  (4.40) 


and  2  '|/+=  20. 

Since  'l'(0)~'l'+  as  ^  »  the  boundary  layer  jump  of  grows  at  ,  the 

jump  of  the  layer  of  n  tends  to  D(0)  and  the  jump  of  the  layer  of  p  tends 
to  zero  exponentially  (see  (4.39)).  The  width  of  the  depletion  layer  which  is  given 
by  (4.29)  grows  logarithmically  with  +  for  constant  ~K.  Figures  17,18,19,20  show 
the  growth  of  the  layer  of  UT¥.  Figure  21  shows  the  boundary  layer  solution  n  for 
U=  -19.194V  (reverse  bias).  The  depletion  region  is  almost  free  of  carriers 
(since  n~l  +  n,  p  ~  0)  . 


0.01 
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0.00  0.12 
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By  proceeding  as  in  Fife  (1973)  (proof  of  leitma  2.1)  we  find  that  (A4),(A5)  has 

a  unique  monotone  solution  V  =  V .  fulfilling  (for  >  C+  >  0) 

¥+ 

0  <  V .  (y)  <  exp((-l+6)y)  ,  y  >  const.^n  ♦  (A6) 

¥+  -  - 

where  0  <  6  <  1  and  the  constant  depends  on  6  only.  The  linearization  of 

2  2 

?(&+>  T',  ’)  (viewed  as  operator  from  C  ([0,1])  into  C([0,l])xK  )  at 
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(x,X)=V^  (^)  is  given  by 


>2 . 


CL  (V.  ,  X)y)(x)  =  (X  y"(x)-g,  (^+  ,  V,  (x,X))y(x)  ,  y(0)  ,  y(l)) 

y+  ¥+ 

In  Ci([0,l])  we  take  ||u||  =  EX*||u^| 


(A7) 


j=0 

to  Fife  (1973),  Lemma  2.5,  we  get 


40,1] 


for  the  following.  Analogously 


vjr 

l|L  lcr  ,  X)  ||  <  const.  X  1e  + 
v+ 


(A8) 


where  the  constant  is  independent  of  >  ,  .  Moreover,  the  linearization  L 


fulfills 


S'i'. 


||L(U,7\)-L(V,A)  II  <  const.  *+e  +||U-V 


(A9) 


for  HU-V^  ||,  ||V-Vt  ||  <|Q|<  1. 


By  inserting  V.  into  (A2)  we  get  using  (A6) 


l|F<VX,V^)||  <  exp(-^i) 


(A10) 


From  Spijker  (1972)  we  conclude  that  (A2)  has  a  locally  unique  rolution  v  if 


evC^)  ■  (V  *+***•  S  <  <» 


(All) 


for  some  p  sufficiently  small.  Then 

||v-V^  (x»X)  ||  <  const.  X  *exp(  ^~^>l)  (A12) 

holds.  (All)  follows  immediately  from  (3.9)  and  (3.10)  follows  from  (A12) 
by  using  (Al). 
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Theorem  3.1  is  a  special  case  of  Theorem  3.2  with  ^ (7. ) .  (3.9)  holds, 

therefore  (3.4), (3.5)  follow  immediately. 

Proof  of  Theorem  3.3.  Again  we  substitute 

¥  =  ++(l-v)  (A13) 

and  get  from  (3.15), (3.8) 

V  -  V 

A2v"  =  ^-(l+p(l,?\)e  +  -n(l,A)e  +  )  ,  0  <  x  <  1  (A14) 

v(0)  =  1  ,  v  (1 )  =  0  .  (A15) 


The  reduced  equation  (A14)  has  the  solution  v  =  0.  From  Fife  (1973)  we  conclude 
that  the  semi-reduced  problem 


-  4?  w 

A2w"  =  r-a  -  e  +  )  ,  0  <  X  <  1 

Y+ 

w(0)  =  1  ,  w(l)  -  0 


(A16) 

(A17) 


which  is  obtained  by  setting  A  =  0  in  the  right  hand  side  of  (A14),  has  a 
solution  fulfilling 

w(x,A)  =  0  (exp  (  -^-(-l  +  6)x))  +0(exp(  -^))  ,  c>0,  0  <  &  <  1  (A18) 

if  A  is  sufficiently  small  (depending  on  »►  ).  From  Ringhofer  (1980), (1981)  we 
conclude  that  the  linearization  of  (A16), (A17)  at  w(x,A)  is  uniformly  stable  with 
respect  to  A.  Since  p (1,A)  = 0 (y^A^)  holds  we  get  (3.16)  by  a  simple  perturbation 
argument  (similar  to  the  proof  of  Theorem  3.2). 

Remark .  Theorem  3.3  is  not  formulated  uniformly  in  'i'  . 


For  the  numerical  solution  of  two-point  boundary  value  problems  we  used 
polynomial  collocation  at  Gaussian  points  (see  de  Boor  and  Swartz(1973),  Russell 
(1974)  and  Weiss  (1974)).  This  method  translates  into  a  system  of  (nonlinear) 


difference  equations  (if  the  two-point  boundary  value  problem  is  nonlinear) 
which  has  to  be  solved  numerically  by  an  iterative  method  (for  example  Newton's 
method)  starting  from  an  initial  guess  for  the  solution  (which  has  to  be 
sufficiently  accurate  for  the  iteration  to  converge) . 

The  order  of  convergence  depends  on  the  number  k  of  collocation  points  in  each 

2k 

subinterval.  For  Gauss  points  the  global  error  at  the  grid  points  is  an  0(h  ) 

where  h  is  the  maximal  mesh  size  assuming  that  the  linearized  system  of  difference 
equations  has  a  uniformly  (as  h>0)  bounded  inverse. 

We  used  the  code  COLSYS  (see  Ascher,  Christiansen  and  Russell  (1973)  and 
Ascher  (1981))  as  collocation  solver  with  k=3  or  k  =  4. 

A  problem  affecting  stability  and  consistency  of  the  collocation  scheme  is 
the  choice  of  the  mesh.  COLSYS  is  a  completely  adaptive  code  refining  an  initial 
mesh  automatically  in  order  to  achieve  a  prescribed  error  tolerance.  However,  we 
did  not  use  this  option  for  solving  boundary  layer  problems.  Instead  we  prescribed 
a  mesh  which  is  exponentially  graded  in  the  layer  region  and  coarse  outside  (where 
the  solution  is  close  to  reduced  solution)  (see  Ascher  and  Weiss  (1981)).  This 
mesh  equidis tributes  the  local  discretization  error  and  employs  only  a  reasonable 
number  of  mesh  points  (error  estimates  and  estimates  on  the  number  of  mesh  points 
are  given  in  Markowich  and  Ringhofer  (1981)  for  a  similar  class  of  problems). 

The  boundary  layer  equations  (4.25), (4.26)  were  solved  by  cutting  the  infinite 
interval  at  a  finite,  far  out  point  and  by  substituting  an  additional  boundary 
condition  at  this  point  for  the  condition  at  infinity  (see  Markowich  (1981)  and 
Markowich  and  Ringhofer  (1981)). 
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